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The understanding of strong-field dynamics near black-hole horizons is a long-standing and challenging prob- 
lem in general relativity. Recent advances in numerical relativity and in the geometric characterization of black- 
hole horizons open new avenues into the problem. In this first paper in a series of two, we focus on the analysis 
of the recoil occurring in the merger of binary black holes, extending the analysis initiated in [H with Robinson- 
Trautman spacetimes. More specifically, we probe spacetime dynamics through the correlation of quantities 
defined at the black-hole horizon and at null infinity. The geometry of these hypersurfaces responds to bulk 
gravitational fields acting as test screens in a scattering perspective of spacetime dynamics. Within a 3 + 1 
approach we build an effective-curvature vector from the intrinsic geometry of dynamical-horizon sections and 
correlate its evolution with the flux of Bondi linear momentum at large distances. We employ this setup to study 
numerically the head-on collision of nonspinning black holes and demonstrate its validity to track the qualita- 
tive aspects of recoil dynamics at infinity. We also make contact with the suggestion that the antikick can be 
described in terms of a "slowness parameter" and how this can be computed from the local properties of the 
horizon. In a companion paper 0] we will further elaborate on the geometric aspects of this approach and on 
its relation with other approaches to characterize dynamical properties of black-hole horizons. 

PACS numbers: 04.30.Db, 04.25.dg, 04.70.Bw, 97.60.Lf 



I. INTRODUCTION 

Understanding the dynamics of colliding black holes (BHs) 
is of major importance. Not only is this process one of the 
main sources of gravitational waves (GWs), but it is also 
responsible for the final recoil velocity (i.e., "kick") of the 
merged object, which could play an important role in the 
growth of supermassive BHs via mergers of galaxies and 
on the number of galaxies containing BHs. The recoil of 
BHs due to anisotropic emission of GW has been known for 
decades Q H] and first estimates for the velocity have been 
made using approximated and semianalytical methods such as 
aparticle approximation @-0], post-Newtonian methods |0- 
and the close-limit approximation lfl2i [l3ll . However, it 
is only thanks to the recent progress in numerical relativity 
that accurate values for the recoil velocity have been com- 
puted GHH- 

Indeed, simulations of BHs inspiralling on quasicircular 
orbits have shown, for instance, that asymmetries in the 
mass can lead to kick velocities < 175km/s ifTii [l5ll . 
while asymmetries in the spins can lead respectively to 
«k < 450km/s or < 4000 km/s if the spins are 
aligned lfl7l [l8l l2(i[| or perpendicular to the orbital angular 
momentum l23h (see l24T.l25ll for recent reviews). 

In addition to a net recoil, many of the simulations show an 
"antikick," namely, one (or more) decelerations experienced 
by the recoiling BH at late times. In the case of merging BHs, 
such antikicks seem to take place after a single apparent hori- 
zon (AH) has been found 12611 (see Fig. 8 of Ref. 112011 for some 
examples). An active literature has been developed over the 
last few years in the attempt to provide useful interpretations 
to this process l26l - [30ll . Interestingly, some of these works do 
not even require the merger of the BHs. As pointed out in If3l"ll 
when studying the scattering of BHs, in fact, the presence of 
the common AH is not a necessary condition for the antikick 



to occur. Furthermore, as highlighted in lL32tl . it is also possi- 
ble to describe this process without ever discussing BHs and 
just using the mathematical properties of the evolution of a 
damped oscillating signaQ. 

Although the presence of a common AH is not a neces- 
sary condition for the appearance of an antikick (which could 
indeed be produced also by the scattering of a system involv- 
ing one or two neutron stars), when a common AH is present 
through the merger of BH binary, we can use information on 
the latter togain insight in the physical mechanisms behind 
the antikiclo We believe that constructing an intuitive picture 
of the dynamics of general relativity in a region of very strong 
field is not only of general interest but also of practical use 
to explain this process. In Ref. GL m fact, a new conjecture 
was suggested in which the antikick produced in the head-on 
collision of two BHs with unequal masses was understood in 
terms of the dissipation of the AH intrinsic deformation. As 
shown in the schematic cartoon in Fig.[T](cf. Fig. 1 of Gil and 
also Fig. [TT] for a comparison with numerical data), the kick 
and antikick can be easily interpreted in terms of simple dy- 
namical concepts. Initially the smaller BH moves faster and 
linear momentum is radiated mostly downwards, thus leading 
to an upwards recoil of the system [stage (1)]. When a sin- 
gle AH is formed at the merger, the curvature is higher in the 
upper hemisphere of the distorted BH and linear momentum 
is radiated mostly upwards leading to the antikick [stage (2)]. 
The BH decelerates till a uniform curvature is restored on the 
AH [stage (3)]. The qualitative picture shown in the cartoon 



On the other hand, if an exponentially damped oscillating signal is present, 
this is indeed a signature of the presence of a BH ringing down. 
Here and in the companion paper we will show that even when a horizon 
is not present, the considerations made here can be extended on a suitably 
defined 2-surface. 
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FIG. 1. Generation of the antikick in the head-on collision of two 
unequal-mass Schwarzschild BHs as described in Qj]. Initially the 
smaller BH moves faster and linear momentum is radiated mostly 
downwards, thus leading to an upwards recoil of the system [stage 
(1)]. At the merger the curvature is higher in the upper hemisphere 
of the distorted BH (cf. red-blue shading) and linear momentum is 
radiated mostly upwards leading to the antikick [stage (2)]. The BH 
decelerates till a uniform curvature is restored on the AH [stage (3)]. 
This cartoon should be contrasted with the results of numerical sim- 
ulations in Fig, II II 



was then investigated by exploiting the analogy between this 
process and the evolution of Robinson-Trautman (RT) space- 
times i33l l34ll and by showing that a one-to-one correlation 
could be found between the properties of the AH perturbation 
and the size of the recoil velocity Jl|] . 

In this paper and in its companion paper (hereafter paper 
I and II, respectively), we provide further support to the con- 
jecture proposed in ^ by extending our considerations in [ 1] 
to more generic initial data and, more importantly, by investi- 
gating in detail numerical spacetimes describing the head-on 
collision of two BHs with unequal masses. To do this we in- 
troduce a cross-correlation picture in which the dynamics of 
the spacetime can be read off from two "screens" provided 
naturally by the black-hole event horizon £ and by future null 
infinity In practice, using the standard 3 + 1 approach 
in general relativity, we replace these screens with effective 
ones represented, respectively, by a dynamical horizon H + 
and by a time like tube B at large spatial distances. We then 
define a phenomenological curvature vector Kf s (t) in terms 
of the (mass multipoles of the) Ricci scalar curvature 2 R at 
H + and show that this is closely correlated with a geometric 
quantity {dPf /dt)(t), representing the variation of the Bondi 
linear momentum in time on J? + . This construction, which 
is free of fitting coefficients and valid beyond the axisymmet- 
ric scenario considered here, correlates quantities on the AH 
with quantities at large distance, thus providing us with two 
important tools. Besides confirming the association of recoil 
dynamics with the dissipation of anisotropic distribution of 
curvature on the AH, it opens a new route to the analysis of 
strong-field effects in terms of purely local quantities evalu- 
ated either on the AH or on other suitable surfaces. 

This first article is organized as follows. Section [Tj] in- 



troduces an executive summary, where the main concepts 
presented in these two papers are summarized for those not 
wishing to enter into the mathematical details. Section [Hi] 
on the other hand, extends the analysis carried out in JT1 
for RT spacetimes by considering more general initial data 
and by analyzing aspects of the evolution of the AH cur- 
vature. Section HVl extends the methodology and diagnostic 
tools to BH spacetimes representing the head-on collision of 
unequal-mass BHs. In particular, we develop the mathemat- 
ical tools necessary to measure the relevant quantities on the 
two screens and we show how they closely correlate. Finally, 
the conclusions are discussed in Sec. [V] while the Appendix 
is used to provide details on our definitions of the correlation 
and matching of time series. 

This paper also builds on the material presented in its com- 
panion paper II, where we present a more detailed discussion 
of the mathematical aspects of our framework. In particular, 
we revisit there the evolution of relevant geometric objects on 
the AH and introduce preferred null normals on a dynamical 
horizon. In our discussion of a newslike function on the dy- 
namical horizon and its relation to the problem of quasilocal 
linear momentum, we also stress the importance of the inner 
horizon when evaluating fluxes across the horizon. 

We use a spacetime signature (— , +, +, +), with abstract 
index notation (first letters, a, b, a.., in Latin alphabet) and 
Latin midalphabet indices, k..., when making explicit the 
spacelike character of a tensor. Greek indices are used for ex- 
pressions in particular coordinates systems. We also employ 
the standard convention for the summation over repeated in- 
dices. Finally, all the quantities are expressed in a system of 
units in which c = G = M Q = 1, unless otherwise stated. 



II. THE CROSS-CORRELATION APPROACH: AN 
EXECUTIVE SUMMARY 

As mentioned above, this section is meant to provide a 
general summary of the results and methodology of papers I 
and II, focusing mostly on the conceptual aspects and leaving 
aside the mathematical details, which can instead be found in 
the corresponding main texts. 

We start by recalling that Ref. ^ suggested an approach to 
study the near-horizon nonlinear dynamics of the gravitational 
fields based on the systematic analysis of the deformations in 
the BH horizon geometry. In particular, it was shown how the 
gravitational dynamics responsible for the antikick after a bi- 
nary merger can be understood in terms of the anisotropies in 
the intrinsic curvature of the AH of the resulting merged BH. 
Considering a RT spacetime, the kick velocity constructed 
from the Bondi momentum (a geometric quantity at null infin- 
ity) was put in a one-to-one correspondence with a quasilocal 
geometric quantity constructed on the horizon, namely, with 
the effective curvature parameter K e g. This geometric pa- 
rameter K e g encodes the part of the AH geometry whose dis- 
sipation through gravitational radiation can be related to the 
final value of the kick. Stated differently, very different bi- 
nary systems, e.g., with very different mass ratio, give rise to 
the same final kick velocity as long as they share the same 
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value of the K c ff parameter. 

The following criteria were employed in Ref. J2l for the 
construction of the curvature parameter K c g: i) K c s should 
not depend on how the AH is embedded in the spacetime; ii) 
K e g should change sign (i.e., it should be an odd function) 
under reflection with respect to a plane normal to a given axis. 
From the first requirement, K e g was constructed in terms of 
the intrinsic geometry of the AH, namely as a functional on 
the Ricci scalar 2 R associated with the induced metric on the 
AH. The ansatz for K c g in Ref. |Q]], compatible with require- 
ment ii) above and within axisymmetry, had the following 
structure 



Lcff 



/oven {M 2 e) x / odd (M 2 e+i) 



(1) 



where Mg's are the so-called isolated-horizon mass multipoles 
associated with a spherical harmonic decomposition of 2 R in 
the axisymmetric case lf35L f36h - The odd part / Q dd accounts 
for the directionality of the kick, whereas the even part / even 
controls its intensity. 

In order to validate this suggestion, we analyzed a family 
of Robinson-Trautman (RT) spacetimes l33l[34ll . representing 
an (eternal) BH together with purely outgoing gravitational 
radiation. The mathematical properties of this class of exact 
solutions is already well understood ll37l - l40ll . therefore this 
spacetime is a good test for numerical schemes i4ll l42ll and 
it is an excellent toy model for problems dealing with radi- 
ation in BH environments 1431 - 15411 . Although the associated 
BH horizon is stationary, these RT spacetimes also contain a 



white-hole horizon H 
outer-trapping horizon 



4( 1 15 51 - 15 7H . or, more precisely, a past 
5811 . whose dynamics offers a particu- 
larly well-suited scenario to test our geometric approach. This 
is shown in Fig. [2] which reports a Carter-Penrose diagram for 
the RT spacetime (see also ll40ll57ll ). The solutions exist for 
u ^ uq and the white hole emits GWs until the Schwarzschild 
spacetime is achieved asymptotically. In practice, numerical 
simulations run up to a finite w nna i and show the exponential 
convergence to a solution which is essentially stationary. 

In Ref. Ql, the functions / even and / Q dd appearing in 
Eq. (Q]) were written in the simplest possible form, i.e., as a 
linear expansion in M/s 

K cB = (a 2 M 2 + a 4 M 4 + . . .) x (a 3 M 3 + a 5 M 5 + ...). 

(2) 

Then, using suitably defined initial data, a set of numerically 
fitted coefficients aj's was found so that a one-to-one depen- 
dence between the final kick velocities «k and K c s at a given 
retarded time u could be found: i.e., Avk = — v(u) = 
A x K c g(u), where v^(u) is the recoil velocity at time u and 
A is a constant. This injectiveQ relation between K c g and 
Vk permits us to understand the degeneracy of the latter, as 
a function of the mass ratio in terms of AH quantities at a 
given (initial) time u (cf. Fig. 3 in Ref. IH and Fig.|9]below). 




3 Note that the relation is not only injective, but also linear. This is ultimately 
due to the writing of -fC c ff as the product of two functions (of even and odd 
multipoles), such that each of these functions is linear in the multipoles. 



FIG. 2. Carter-Penrose diagram for the RT spacetime 14(115711 . The 
solutions exist for u ^ Uo. The white hole emits GWs until the 
Schwarzschild spacetime is achieved, although the numerical simu- 
lations run until a finite Ufinai, when an essentially stationary solution 
is found. 



Moreover, the good quantitative agreement between i>k cal- 
culated from full binary BH numerical simulations and from 
RT models, suggested the presence of a generic behavior in 
this physical process. Overall, therefore, the work in Ref. |Q]] 
provided an approach to understand global recoil properties 
in terms of (quasi-)local quantities on the AH, and an intuitive 
guideline to interpret the black-hole recoil properties in terms 
of the dissipation of AH geometric quantities. 

Despite the valuable insight, the treatment presented in 
Ref. HI] had obvious limitations. First, the ansatz for K e g 
in Eq. (Q} is not straightforwardly generalizable to the nonax- 
isymmetric case. Second, the phenomenological coefficients 
a/s in Eq. (fJJ depend on the details of the employed RT initial 
data. Finally, the white-hole horizon analysis in RT space- 
times needs to be extended to the genuine BH horizon case. 
All of these restrictions are overcome in the work reported in 
papers I and II. 

While the focus in Ref. 0]] was on expressing the difference 
between the final kick velocity and the instantaneous kick 
velocity v^u) at an (initial) given time u, in terms of the ge- 
ometry of the common AH at that time u, we here focus on ge- 
ometric quantities that are evaluated at a given time during the 
evolution. More specifically, we will consider the variation 
of the Bondi linear momentum vector in time (dP® / du)(u) 
as the relevant geometric quantity to monitor at null infinity 
J^ + . To this scope, we need first to construct a vector K l oS (v) 
(function of an advanced time v) as a counterpart on the BH 
horizon H + . Then, we need to determine how K l cS (v) on H + 
correlates to {dP® / du)(u) at J- + . 

In the RT case, the causal relation between the white-hole 
horizon H~ and null infinity made possible to establish 
an explicit functional relation between dv^/du and K' cS (u). 
In the case of generic BH horizon, however, such a direct 
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FIG. 3. Carter-Penrose diagram illustrating the scattering approach 
to near-horizon gravitational dynamics in a generic spherically sym- 
metric collapse. The event horizon $ and null infinity pro- 
vide spacetime canonical screens on which geometric quantities, re- 
spectively accounting for horizon deformations and wave emission, 
are defined. Their cross-correlation encodes nontrivially information 
about the bulk spacetime dynamics. 

causal relation between the inner horizon and is lost 
(compare Fig. [2] with Figs. [3] and |4j. However, since their 
respective causal pasts partially coincide, nontrivial correla- 
tions are still possible and expected. This can be measured 
through the cross-correlations of geometric quantities h[ nn (v) 
at % + and h out (u) at J? + , both considered here as two time 
seriefl In particular, we will take K l eS (v) as h- mn (v) and 
(dPf/du){u) as h out (u). 

This approach to the exploration of near-horizon gravita- 
tional dynamics resembles therefore the methodology adopted 
in scattering experiments. Gravitational dynamics in a given 
spacetime region affects the geometry of appropriately chosen 
outer and inner hypersurfaces of the BH spacetime. These 
hypersurfaces are then understood as test screens on which 
suitable geometric quantities must be constructed. The corre- 
lations between the two encode geometric information about 
the dynamics in the bulk, providing information useful for 
an inverse-scattering approach to the near-horizon dynamics. 
As a result, in asymptotically fiat BH spacetimes, null in- 
finity and the (event) BH horizon S 1 provide preferred 
choices for the outer and inner screens. This is nicely sum- 



Note that the meaningful definition of time series cross-correlations re- 
quires the introduction of a (gauge-dependent) relation between advanced 
and retarded time coordinates v and u. In an initial value problem this is 
naturally provided by the 3 + 1 spacetime slicing by time t. 



FIG. 4. Carter-Penrose diagram for the scattering picture in a Cauchy 
initial value approach. The dynamical horizon H + and a large- 
distance timelike hypersurface B provide inner and outer screens. 
Note that a the dynamical horizon is split in two portions: outer and 
inner (solid and dashed blue lines, respectively) and that the 3 + 1 
slicing sets a common time t for cross-correlations. 

marized in the Carter-Penrose diagram in Fig. [3] which illus- 
trates the cross-correlation approach to near-horizon gravita- 
tional dynamics. The event horizon <§ and null infinity 
provide natural spacetime screens on which geometric quan- 
tities, respectively, accounting for horizon deformations and 
wave emission, are defined. Their cross-correlation encodes 
information about the bulk spacetime dynamics. 

Although the picture offered by Fig. [3] is quite simple and 
convincing, it is not well adapted to the 3+1 approach usually 
adopted in numerical studies of dynamical spacetimes. In- 
deed, neither the BH event horizon nor null infinity are in gen- 
eral available during the evolution^. However, we can adopt 
as inner and outer screens a dynamical horizon H + (future 
outer- trapping horizon 15 8M60IO and a timelike tube B at large 
spatial distances, respectively. In this case, the time function t 
associated with the 3 + 1 spacetime slicing provides a (gauge) 
mapping between the retarded and advanced times u and v, 
so that cross-correlations between geometric quantities at H + 
and B can be calculated as standard time series hi nn (t) and 
h ut{t)- This is summarized in the Carter-Penrose diagram in 
Fig. 2J which is the same as in Fig.[3j but where the 3 + 1 slic- 
ing sets an in-built common time t for cross-correlations be- 
tween the dynamical horizon H + (i.e., the inner screen) and a 
large-distance timelike hypersurface B (i.e., the outer screen). 



5 The latter would properly require either characteristic or a hyperboloidal 
evolution approach. 
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Within this conceptual framework it is then possible to de- 
fine a phenomenological curvature vector Kf (t) in terms 
of the mass multipoles of the Ricci scalar curvature 2 R at 
H + and show that this is closely correlated with a quantity 
(dPf /dt) (t) on B, representing an approximation to the vari- 
ation of the Bondi linear momentum time on How to 
do this in practice for a BH spacetime is the subject of the 
following sections. 



III. ROBINSON-TRAUTMAN SPACETIMES: A TOY 
MODEL 

We recall that the RT spacetimes are a class of solutions 
of the vacuum Einstein equations admitting a congruence of 
null geodesies which are hypersurface-orthogonal, shear-free 
but with nonvanishing expansion. As such, it can be regarded 
as a white hole emitting GWs, thus representing a valuable 
tool for studying the spacetime geometry in physical condi- 
tions that are similar to the final stages on the dynamics of BH 
binaries 13311 . The RT metric can be written as 14911 



ds 2 



K 



2M n 



2rd u Q 

Q 



du — 2dudr 



-dQf 



(3) 

where Q = Q(u, ft), u is a null coordinate, r is an affine pa- 
rameter of the outgoing null geodesies, and dil 2 = d9 2 + 
sin 2 9 dip 2 is the metric of a unit sphere S 2 . Here 
is a constant and is related to the mass of the asymptotic 
Schwarzschild BH, while the function K(u, f2) is the Gaus- 
sian curvature of the surface with r = 1 and u = const., and 
is given by 



K(u,Q) = Q 2 (l + A 2 J lnQ), 



(4) 



being the Laplacian operator on the unit sphere S 2 . The 
Einstein equations then lead to the RT evolution equation 



d u Q(u, fi) = -Q- 



12M X 



(5) 



Any regular initial data Q = Q(uo, f2) will smoothly evolve 
according to (f5J) until it reaches a stationary configuration cor- 
responding to a Schwarzschild BH at rest or moving with a 
constant speed l39ll . Equation © implies the existence of the 
constant of motion A = § s dQ/Q 2 , which clearly represents 
the area of the surface u = const., r = 1 and can be used to 
normalize Q so that A = 4ir. 

The dynamical compact object modeled by RT spacetimes 
is described by the past AH, which has a vanishing expansion 
of the ingoing future-directed null geodesies. Such past AH is 
described by the surface r — R(u, f2) satisfying Il55l - D7ll 



An lnR = K - 



2M a 
R 



(6) 



The line element restricted to the AH surface u — const, 
and r = R{VL) is 



(7) 



which, from Eq. (|6), has a Gaussian curvature 

2M na 



A' = 



R? 



(8) 



On the other hand, the mass and momentum are computed 
at future null infinity using the Bondi 4-momentum as ll3~7ll38[ 



Moo 
47T 



Q~ 3 



dQ,. 



(9) 



with {77"} = {1, sin cos 93, sinflsin^, cos#}. 

From now on, we restrict our problem to axisymmetry and 
introduce x = cos 9. Clearly, all the physically relevant in- 
formation is contained in the function Q(u,x), and this in- 
cludes also the gravitational radiation, which can be extracted 
through the radiative part of the Riemann tensor l33l l34ll . 
which in axisymmetry is given by 



r* 4 



On 



(l-x 2 )d 2 Q 



Q 



(10) 



The dynamics of this solution can be summarized in Fig. [2] 
which shows the Carter-Penrose diagram for the RT space- 
time. The final configuration is a stationary nonradiative solu- 
tion which has the form 14911 



S±{9) = Q(oo, 0) 



(1 =F Vkx) 



(11) 



Note that since the Bondi 4-momentum of the stationary solu- 
tion is 



P a (oo) 



{i,o,o,±« k }, 



(12) 



the parameter in Eq. (fTTT > is interpreted as the velocity of 
the Schwarzschild BH in the z-direction. 



A. Mass multipoles 

Given a closed 2-surface S, the invariant content of its in- 
trinsic geometry is encoded in the Ricci scalar curvature 2 R 
associated with the induced metric q a b on S. Moreover, if S 
is an axisymmetric surface, with <\> a as the axial Killing vector, 
a preferred coordinate system (9, (p) can be constructed such 
that q a b has the form ll35l, l36tl 

q ab dx a dx b = R 2 H (/^sin 2 ^ d9 2 + f d((A , (13) 

where f(9) = q ab (j) a cj) b / R 2 , with R H the areal radius (A = 
J s dA = 47ri? 2 ). The coordinate 9 is determined by 



(14) 



where the coordinate £ is defined by £ = cos6> and \ ba is the 
alternating symbol. In addition, the normalization condition 
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(dA = must be imposed. We note that the Ricci scalar 
2 R on S can be written as l35ll 

and that regularity conditions on the metric impose 

Urn /(C) = 0, Urn %(() = ±2 . (16) 

A crucial feature of this coordinate system is that the associ- 
ated expression for the area element is proportional to that 
of the "round sphere" metric dA = R^smB d9d(p. This 
provides the appropriate measure on S to define the standard 
spherical harmonics Yg, TO =o(#) with the standard orthonormal 
relations 



Y efi 0)Y e , iO 0)dA = R 2 Jw , 



(17) 



so that dimensionless geometric multipoles Ig can be intro- 
duced as the spherical harmonics components of the Ricci 
scalar curvature 2 R IB5I1 



1 



2 RY lfi {6) dA, 2 R 



— Y,hYi fi {6). (18) 

H „ =0 



The mass multipoles M/s are then defined as appropriate di- 
mensionful rescalings of the geometric Ig's 



4tt Af H (i? H r 



2n 



2tt 



(19) 



where M H denotes some appropriate quasilocal mass for the 
surface S. Because we will consider here initial data with 
zero angular momentum, M H will denote the irreducible mass 
Mi rr = tJA/JWtt) = Rb_/2. For later convenience, we in- 
troduce the rescaled geometric multipoles Ig 

1 . 4 _ 



Mr 



(20) 



with dimensions [Ig] = [length] 2 . 
ture can then be written as 



The Ricci scalar curva- 



2 R = ^m 



to ■ (21) 



1=0 



A crucial remark for the discussion in Sec. lIV Bl is the vanish- 
ing of the I = 1 mode, i.e., I\ = 0, which can be interpreted 
as a choice of center of mass frame of the AH in lf35ll . This 
follows by first inserting expression <TT3T > into the definition of 
I\, so that 1 1 ex. J_ 1 /"(C) C^C an d then by making use of 
regularity conditions ( TTST ) after integrating by parts. 

In the particular case of a RT spacetime, the preferred ax- 
isymmetry coordinate system (9, (p) is related to the RT spher- 
ical coordinates as (p = cp and 9 = 6(9) satisfying 



smdR{0) 2 

Jr^qW 



(22) 



This equation is solved with the condition C(0) = 1 and then 
one computes the mass multipoles moments through the ex- 
pression 



M e = M 



£+1 



\0)R(6) 



dCl. 



(23) 



B. The numerical setup 



As discussed in detail in Ref. 14911 , for the numerical solu- 
tion of the Einstein equations we introduce a Galerkin decom- 
position for Q(u, x) 



N 

Q(u, x) = ^ bg(u)Pg(x), 

£=0 



(24) 



where Pg(x) stands for the Legendre polynomial of order I. 
By using standard projection techniques, Eq. (O can be writ- 
ten as a system of ordinary differential equations 



2£ + 1 



(Q 3 x [(l-x 2 )d x K] ,Pg] 



t = 0,l,.-.,N, 
(25) 

where the inner product is given by (/, g) = j ' fgdx. In 
this way, the Cauchy problem for the RT Eq. ( T25l l consists 
basically in choosing the initial value of the mode functions 
bg(u) according to 



6/(0) = 



2£+ 1 



(Q(0,x),P/>, 



(26) 



and then to solve the initial value problem given by d25l >. Note 
that, asu-4 oo, bg — > for £ > 1 and that the nonzero modes 
must satisfy fei(oo) 2 — &o(°o) 2 = 1> with tne nna l v k param- 
eter of Eq. dTTT > being given simply by = — &i(oo) /6o(oo). 

Equation (O can be solved for the horizon R(u, 9) either 
by imposing regular conditions on the boundary dgR(0) = 
and dgR(n) = and using an ordinary shooting method 
to find i?(0) and R(ir), or by following the approach in ll6lll 
introducing another Galerkin decomposition on the horizon 
and truncated at the order Na 



In 



R 



2M„ 



= y]cj(u)Pe(x). 



(27) 



The projection of Eq. (|6]l on the basis of the Legendre polyno- 
mials couples the known Galerkin modes bg with the unknown 
coefficients eg and the resulting algebraic nonlinear system 
can be easily solved via a Newton-Raphson method. 



C. The initial data 

In general, any family of regular functions Q, i.e., < 
Q(uq, x) < co, Vx G [— 1, 1] can be used as an initial data for 
the RT spacetime. For any of such family, one can set a param- 
eter Qq to ensure the constant of motion to be A = 47r. More- 
over, the deformed BH will not be initially at rest in general. 
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FIG. 5. AH of the head-on initial data given in Eq. {28}. The parameters q and w control the deformation of the surfaces. The final configuration 
is a Schwarzschild BH (continuous black line) and it is achieved after the deformation is dissipated with the emission of GWs. Our results 
confirm the interpretation given li48h for q as the mass ratio of the BHs. However, we attach no physical meaning to w as done in [48], for RT 
models a single deformed horizon, and nothing can be said about the velocities of the individual colliding BHs. 



As a result, given the initial velocity «k.o = -P 3 (0)/P°(0), 
we perform a boost P° = A Q J g(i?k,o)-f >/3 with A a p(vk,o) the 

3 

associated Lorentz transformation, so that P (0) = by con- 
struction. The kick velocity is then defined at any time as 
v{u) = P 3 (u)/T^(u). 

Despite this overall simplicity, the definition of initial data 
that is physically meaningful represents one of the main diffi- 
culties in the study of RT spacetimes. Since we are here more 
interested in a proof of principle than in describing a realis- 
tic configuration, we have adopted both a prescription remi- 
niscent of a "head-on" collision of two BHs 14811 and a new 
variant of it. 



1. "Head-on" initial data 

As a first set of initial data we consider the one developed 
in Ref. H 



Q(0,0) = Qa 



1 



<1 



1 -2 



\f\ — WX + wx 



(28) 



which was interpreted to represent the final stages (i.e., af- 
ter a common AH is formed) of a head-on collision of two 
boosted BHs wit h op posite velocities w € [—1, 1] and mass 
ratio q <G [0,1] 14811 . Figure [5] shows the shape of the sur- 
face R(uq, 9) for different values of those parameters. It is 
worth remarking that despite the name, this initial data does 
not represent a binary system but is, strictly speaking, only a 
distorted horizon. In a more conservative approach, one can 
regard q and w just as free parameters that control the de- 
formation on the horizon, and this is the view we will adopt 
hereafter. However, a number of interesting analogies with the 
head-on collision of two BHs have been suggested ll48[ 
and will be further discussed below. 



The interpretation of q as a mass-ratio parameter is not to- 
tally unreasonable. For instance, Fig. 3 of [1] showed the final 
value of the velocity in RT spacetime evolved from the head- 
on initial data against the reduced mass ratio v = q/(l + q) 2 . 
The curve obeys the distribution 



Av 2 sJ\ -±v(l + Bv) 



(29) 




as found in all numerical simulations H24H . with the value of 
A and B depending on the particular choice of w. 

Since the solution does not exist for u < uo, it is impossi- 
ble to assign a value for w that could account for any previous 
stages on the evolution of the binary. With the original in- 
terpretation of w as the velocities of the BHs, a first trivial 
estimate, as proposed by lO . is to assume a Newtonian evo- 
lution of two particles with masses M\ and Mi, which start 
at rest at an initial distance Lq. At a given distance L, in the 
frame where v\ = — v-i = w one has 



(30) 



with M = Mi + M 2 . Choosing L ~ 6 M and L ~ 2 M, one 
obtains w ~ 0.41. Furthermore, still in Ref. |Q]] it was shown 
that w = 0.425 presents a surprisingly good match with some 
results found in the close-limit approximation, where the ini- 
tial data for the ringdown phase was given by a previous 
plunge with the BHs inspiralling towards each other from the 
innermost circular orbit until ~ 2M M26I1 . 

It is important to remark that one should not expect a com- 
plete agreement between the values of A and B from the 
head-on collision in RT spacetimes with the ones found in 
numerical-relativity calculations of binary BHs in quasicir- 
cular orbits lfl5ll . The first ones, in fact, (and modulo the 
interpretative issues discussed above) can only account for 
the post-merger phase, while the second ones account for the 
whole recoil. A complete discussion on the dependence of A 
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FIG. 6. Left panel: Horizon mean curvature K H for w = 1/2 shown as a function of reduced mass ratio v and of the polar angle x = cos 6. 
Note that the distribution is symmetric with respect to the equatorial plane for v = and v = 1/4. i?/g/if panel: Horizon curvature for some 
representative values of the reduced mass ratio. The low-;/ branch is characterized by large curvature gradients across the AH but small values 
of the curvature, while the high-;/ branch small curvature gradients and large values of the curvature y]]. 



and B with respect to w can be found in ll52ll . For the sake of 
convenience we will fix w = 0.5, but our results do not de- 
pend upon this choice. The substitution w — s- —w just changes 
the sign of the recoil velocity. 

The relation between the kick velocity and the reduced 
mass ratio expressed by Eq. d29b has a peak ug 13 * for v max sa 
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FIG. 7. Lower-order mass moments for w — 1/2. Note that they 
all vanish at u — 0, while only the odd ones are zero for v = 1/4. 
As in Fig. 3 of jH], the colors represent the two different branches 
composing the curve versus v. The color changes at the mass 
ratio for which the recoil velocity is at a maximum. 



0.195 III5I1 , As the recoil vanishes for v = and v = 1/4, 
there will always be two values of the mass parameter leading 
to the same recoil when v < u™ ax . It is natural to wonder 
whether two completely different systems (namely systems 
with different mass ratios) share some common physical prop- 
erty which could lead to the same recoil. Thinking in terms of 
the horizon's intrinsic deformation provides a simple way to 
explain this degeneracy. System with v < i/ max , in fact, are 
characterized by large curvature gradients across the AH but 
small values of the curvature, while system with v > j/ max are 
characterized by small curvature gradients and large values of 
the curvature. This intuitive picture can be best appreciated 
in Fig. [6] which reports the horizon mean curvature K H for 
w = 1/2 as a function of reduced mass ratio v and of the 
polar angle x = cos 9 (left panel) or for some representative 
values of the reduced mass ratio (right panel). Note that the 
low-^ branch is characterized by large curvature gradients be- 
tween the north and south poles of the AH, but small values 
of the curvature, while the high-i^ branch is characterized by 
small curvature gradients and large values of the curvature. As 
discussed in Ref. QJJ] , it is the product of the deformations on 
the horizon with the gradients across the equator that yields 
the same recoil for two apparently different systems. 

Complementary information to the one in the left panel of 
Fig.|6]is depicted in Fig. [7] which shows the typical behavior 
of the lower-order mass moments as a function of the reduced 
mass ratio. Notice they all vanish for v = 0, since this con- 
figuration represents an undistorted BH. For v = 1/4, on the 
other hand, only the odd modes are zero, indicating that the 
configuration is symmetric with respect to the equatorial plane 
and the emission of GW will not give rise to a recoil. Also 
note that the maximum of the odd modes does not correspond 
to the mass ratio at which the recoil velocity reaches its high- 
est value. Even though a net emission of momentum will only 
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FIG. 8. Left panel: Curvature evolution for the distorted- AH initial data with w — 1/2, v = 0.11. Note the small curvature excess at the south 
pole, which is rapidly radiated away to yield an almost uniform distribution after u/M^o ~ 1.5. Right Panel: evolution of the recoil velocity 
for representative values of the parameter g, some of which lead to no-monotonic changes in the velocity. 



take place when there is an asymmetry on the horizon across 
the equator (i.e., M dd 7^ 0), the intensity of the emission will 
also depend on how deformed the BH is (i.e., M even ^ 0). 

Without loss of generality, we can use the even modes to 
measure overall distortions on the horizon, while the odd ones 
measure the asymmetries between the north and south hemi- 
spheres. To account for both contributions, we constructed 
in an effective-curvature parameter as the product of two 
functions depending solely on the even or the odd modes, 
i.e., A'off = / cvon (M 2n ) x /odd (M 2 „+i). This quantity rep- 
resents a measure of the global curvature properties of the ini- 
tial data, from which the recoil depends in an injective way. 
Indeed, Fig. 4 in Jl|] showed that with a suitable choice of 
coefficients, i.e., A' ff = M%\ J2n=i -^2n+i/3 n_1 |, the cor- 
relation between K e g measured at the initial time against the 
final velocity is actually linear. 

2. Distorted-AH initial data 

The evolution produced by the head-on initial data (l28l l 
leads to a monotonic increase/decrease of the recoil velocity 
once the initial data is specified on the white hole. However, 
a more complex (i.e., nonmonotonic) dynamics can be eas- 
ily produced through a simple variation of the head-on initial 
data. We refer to this new family of initial conditions as to the 
"distorted AHs" and we express them as 

Q(uq, x) = Qho(x; q, w) + qx 2 Qn (x; q, -w), (31) 



where Qno(x;q,w) corresponds to the head-on initial data 
(l28l l. Clearly, this is just a mathematical choice and no physi- 
cal significance can be associated to this initial data. 

Note that OTT l maintains the symmetries provided by the 
mass-ratio parameter q: for q = one recovers the nonde- 
formed Schwarzschild BH and q = 1 gives an even initial 
data, i.e., Q(uq,x) = Q{uq,—x), which leads to a zero fi- 
nal recoil. Furthermore, note that the resulting recoil veloc- 
ity does not lead to the scaling expressed by Eq. (|29l , thus 
giving strength to the idea that the head-on initial data is 
closely related to the merger of a binary system as proposed 
inRefs. EHH. 

As anticipated, the use of this initial data leads to a more 
interesting dynamics and this is shown in Fig. [8] whose left 
panel reports the curvature evolution for w = 1/2 and v 
0.1 1 . while the right panel reports the evolution of the recoil 
velocity for representative values of the parameter q, some of 
which lead to nonmonotonic changes in the velocity^. Note 
the small curvature excess at the south pole, which is rapidly 
dissipated as GWs are emitted so as to yield an almost uniform 
distribution. As pointed out in |Q]], the velocity reaches its final 
value when there is no asymmetry in the deformation between 
the north (x > 0) and south (x < 0) hemispheres, i.e., after 
u/Moo ~ 0.5. 

To prove that the approach discussed in the previous sub- 
section is indeed generic, we define an effective curvature 
K e ff also for this family of initial data, again in terms of the 
product of odd and even mass moments 



6 This evolution is related to the Bondi momentum as defined by i37l l38ll 
and given by Eq. {9}. Recently, a different approach has been proposed 
in [523, where they showed a slightly different profile for the velocity time 



evolution. However, the difference is not important in our argument, since 
we are mainly concerned with values of the curvature at the initial time, 
when we boost P a to its rest frame, and its correlation with the asymptotic 
final velocity, when the momentum is unambiguously defined. 
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K cff = (M 2 + a 4 Mi + a 6 M 6 + a s M 8 + a w M w ) x (M 3 + a 5 M 5 + a 7 M 7 + a 9 M 9 ) , (32) 
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FIG. 9. Difference Av^(u) = «k(^oo) — «k(w) against the effective- 
curvature parameter for several values of the "initial" time u. The 
linear relation between the effective curvature {32} and the kick ve- 
locity is preserved along the evolution as indicated by the different 
symbols, each representing a specific time in the evolution. This re- 
mains the case also at late times as shown in the inset. 



FIG. 10. Evolution of the velocity (red curve) measured with the flux 
of momentum carried out by the GWs. Note the antikick at about 
t/M « 145 that decelerates the system before the final kick velocity 
is reached. The GW signal is instead shown in the inset, namely the 
dominant $4' multipole (blue curve). 



to find that the set of coefficients 04 = 0.304, a@ = 0.178, 
a 8 = 0.086, 010 = -0.186 and a 5 = 0.076, a 7 = -0.090, 
ag = —0.183 leads to the expected injective, and actually 
linear, behavior. Clearly, and as it is also natural to expect, the 
coefficients are different from those found for the head-on ini- 
tial data, and they will always depend upon the specific family 
of initial data considered. The remarkable feature though, is 
that they remain constant in time. This is illustrated in Fig. [9] 
which shows that the effective curvature ( f32b is still linear 
with respect to the relative velocity Ai> k (w) = w(iioo) — v(u) 
at any time during the evolution (this is shown by the different 
coloured symbols, each of which refers to a specific time) and 
also at late times (see inset). This time-independent property 
is general, and not limited to this particular family of initial 
data. In particular, it is also found, for instance, in the head- 
on case. This result reflects the fact that the deformations of 
the horizon evolve in time in a self-similar manner, so that al- 
though the ranges for K e g change in time (becoming smaller 
as the deformations are radiated away), the corresponding re- 
coil velocities maintain the same proportionality (cf. Fig. [5). 

As a concluding remark, we can summarize as follows the 
insight gained through the study of RT spacetimes: the con- 
struction of the effective-curvature parameter depends quan- 
titatively on the family of initial data considered, but that for 
any choice of data, it is possible to find an explicit expression 



that relates the effective-curvature parameter to the (final) re- 
coil velocity through a one-to-one mapping. What however 
will not depend on the initial data is the functional dependence 
of the effective-curvature parameter as expressed by ([]]). In- 
deed, in the next section we will generalize the idea and func- 
tional form of the effective-curvature parameter to account for 
the dynamics in binary BH spacetimes. 



IV. BLACK-HOLE SPACETIMES: HEAD-ON COLLISIONS 

A. Numerical Setup and Results 

The numerical solution of the Einstein equations has been 
performed using a three-dimensional finite-differencing code 
solving a conformal-traceless "3+1" BSSNOK formulation of 
the Einstein equations (see f? 1 for the full expressions) using 
the Einstein Toolkit 062I1 . the Carpet |63J1 adaptive 
mesh-refinement driver, AHFinderDirect Il64ll to track 
the AHs, and QuasiLocalMeasures i65ll to evaluate the 
mass multipoles associated with them. Recent developments, 
such as the use of 8th-order finite-difference operators or the 
adoption of a multiblock structure to extend the size of the 
wave zone have been recently presented in l66l ? ]. Here, how- 
ever, to limit the computational costs and because a very high 
accuracy in the waveforms is not needed, the multiblock struc- 
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FIG. 11. Realization of the cartoon in Fig. Q] with numerical data from a simulation of head-on collision with mass ratio q — 1/2. The color 
code shows the mean curvature on the apparent horizons, which has the same anisotropic behavior of the intrinsic curvature. As described in 
FigQ] once the common horizon is formed, the curvature is stronger in the region of the smaller BH and is dissipated during the evolution. 



ture was not used. Also, for compactness we will not report 
here the details of the formulation of the Einstein equations 
solved for the form of the gauge conditions adopted. All of 
these aspects are discussed in detail in [? ], to which we refer 
the interested reader. 

Our initial data consists of head-on (i.e., zero angular mo- 
mentum) Brill-Lindquist initial data with a mass-ratio of q = 
1/2. The initial separation of both BHs is 10M and they are 
initially located at (0, 0, 6.6666) and (0, 0, -3.3333) to reflect 
their center-of-mass offset. Both BHs have no angular nor 
linear momentum initially. We use a 3D Cartesian numerical 
grid with seven levels of mesh-refinement for the higher mass 
and eight levels of mesh-refinement for the lower mass BH. 
The resolution of our finest grid is h = M/64, while the an- 
gular grid used to find the AHs and evaluate any property on 
these 2-surfaces has a resolution of 65 points in (^-direction 
and 128 points in ^-direction. The extraction of GWs is per- 
formed calculating ^4 at finite-radius detection spheres with 
radii of r x = 60 M , r 2 = 85 M and r 3 = 145 M and then 
extrapolating to infinity. 

Some of the most salient results of the numerical simula- 
tions are summarized in Fig. [T0l which reports the evolution 
of the recoil velocity (red curve) measured with the flux of 
momentum carried out by the GWs. Note the development of 
the antikick at about t/M w 145 (followed by several smaller 
oscillations) that decelerates the BH before the final kick ve- 
locity is reached. Also shown in the inset is GW signal in 
its larger multipolar component (blue curve). Similarly, 
Fig.[TT]provides a realization of the cartoon in Fig. Q] with nu- 
merical data from a simulation of head-on collision with mass 
ratio q = 1/2. Shown with a color code is the mean curvature 
on the apparent horizons, which shares the same qualitative 
properties, and in particular the anisotropic behavior, of the 
intrinsic curvature. As intuitively described in FigQ] once the 
common horizon is formed, the curvature is stronger in the 
region of the smaller BH and is dissipated as the evolution 
proceeds. Note that the curvature distribution is anisotropic 
already at the beginning, as the BHs are tidally distorting each 
other. 



B. Geometric quantities at the BH horizon: iff (t) 

As discussed in Sees. Q] and H2 the analysis of the recoil dy- 
namics in generic BH spacetimes requires a shift with respect 
to the methodology used in RT spacetimes. When consider- 
ing standard 3 + 1 numerical solutions of BH spacetimes, in 
fact, we study the near-horizon dynamics responsible for the 
BH recoil in terms of the time cross-correlations between a 
vector (dP® / dt)(t) at a large-radius hypersurface B and an 
effective-curvature vector Kf (t) constructed from the intrin- 
sic geometry on the dynamical BH horizon H + . The vector 
(dP® / dt)(t) on B approximates the Bondi linear momentum 
flux (dPP /dt)(t) at J? + . From now on we will systemati- 
cally refer to (dPf /dt){t) (and to instead of B), under- 
standing that we are actually using an approximation. 

The construction of Kf R (i) at H + is based on the following 
two guidelines: a) Kf (t) is built out of the intrinsic geome- 
try Ricci scalar curvature 2 R on H + sections; b) the functional 
form of (dPf / di)(t) in terms of the geometry at guides 
the choice of the functional dependence of Kf (t) on 2 R. The 
first requirement is motivated by the success in the RT case, 
whereas the second one aims at preserving those basic struc- 
tural features of the specific function to be cross-correlated. 

Following these guidelines, we start by expressing the flux 
of Bondi linear momentum at null infinity. In terms of a re- 
tarded time u parameterizing J? + , its Cartesian components 
can be written as 



dP[ 
du 



(u) 



lim 

(u,r— >o 



8tt 



Sl \Af(u,n)\ 2 dn, (33) 



where r parameterizes the large-radius spheres S UtT along a 
u = const, hypersurface, r 2 dfl is the area element on S u ^ r , 
s l is its normal unit vector with Cartesian components s 1 = 
(sin#cos</?, sin(9sin(/3, cos#), and the news functions Af(u) can 
be expressed in terms of the Weyl scalar ^4 as 



Af(u,fi) = / y 4 (u',n)du' . 



(34) 



In our 3 + 1 setting with an outer boundary at a finite spa- 
tial distance we need to express the flux with respect to the 
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time function t parameterizing the spatial slices £ t , so that 
we replace S u . r with St, r 



—r-(t) = lim — — 

at r->oo IQn 



^ A (t',Sl)dt' 



dSl. 



(35) 

and where we can think of t (related to u by u = t — r 
near J? + ) as parameterizing the cuts of by hyperboloidal 
slices or, alternatively, the cuts of the timelike hypersurface B 
approximating J? + at large r. We can now rewrite expres- 
sion (l35T l in terms of a generic vector transverse to S u , r 
(i.e., with a generically nonvanishing component along the 
normal to S u>r ), so that the component of the flux of Bondi 
linear momentum along is 



r^-tt) = lim — — 

(Xt r->oo 167T 



St. 



* 4 (i / ,n)dt / 



(36) 



We take this expression as the starting point for the con- 
struction of Kf . It provides the functional form of the Bondi 
linear momentum flux in terms of the relevant component of 
the Riemann tensor at J? + , namely ^4. Then, the two above- 
mentioned guidelines for the construction of Kf s can be met 
by considering a heuristic substitution of ^4 by 2 R in expres- 
sion (13 6b . 

It is important to note that in the same way in which the 
outgoing null coordinate u parameterizes naturally J^ + , the 
ingoing null coordinate v, which runs along is a natural 
label to parameterize the horizon H + . However, within our 
3 + 1 setting, we use Eq. d36l l as the ansatz leading to the 
following proposal for the component K c ff[£](t) of Kf s (t) 
along a vector (tangent to the slice £f) transverse to the 
section St ofH + 



1 

Ifhr 



St 



dA , (37) 



with 



Af(t, Sl) 



— / 2 



2 R(t',Sl)dt' +J\f t "(Sl). 



(38) 



In the equations above, dA is the area element of St, the 
global negative sign accounts for the relative change of the ori- 
entation of outgoing vector normal to inner and outer bound- 
ary spheres, s l are the components of the unit normal vector 
to St tangent to S t , and J\f tc (Sl) is a generic function on the 
surface to be fixed. 

Some remarks are in order concerning expressions d3Tb and 
d38l . First, there is a clear asymmetry between expressions 
d36*l l and (l37l i when substituting the complex quantity "J 4 at 
J? + (encoding two independent modes corresponding to the 
GW polarizations) by the real quantity 2 R on the inner hori- 
zon (a single dynamical mode). Inspection of Eq. (l36l l im- 
mediately suggests an alternative to 2 R by the natural inner 
boundary analogue of ^4, i.e., ^o- However, this strategy 
must face the issue of identifying an appropriate null tetrad at 
H + for the very construction of ^q. Second, the lower limit 



in the time integration, i ^ -00, appearing in Eq. d36b must 
be replaced by the time t c of first appearance of the common 
horizon, when quantities as 2 R(t, SI) start to be well-defined. 
However, there is still a deeper difference between Af(t, SI) 
and N it, SI). Even though one can construct the former as in 
Eq. d34b , i.e., as the time integral of ^4, the definition of the 
news function is local in time depending only on quantities on 
St and not requiring the knowledge of the past history of 
The latter though is here defined as the time integral of 2 R and 
there is no reason to expect the same local-in-time behavior, 
specially as t — > 00. In particular, we fix the function A/"' c (SI) 
by imposing lim N (t, Sl) = 0. 

t— too 

All the points raised above are addressed in detail in the 
accompanying paper II and we adopt here a purely effective 
approach to K l cS (t), since 2 R represents an unambiguous ge- 
ometric object that captures the (possibly many, if matter is 
included) relevant dynamical degrees of freedom in a single 
effective mode. Ultimately, this heuristic proposal for the ef- 
fective curvature is acceptable only as long as it can be cor- 
related with dPP/dt, and this is what we will show in the 
following. 



1. Axisymmetric BH spacetimes 

As a first application of the ansatz d37b , we consider the 
axisymmetric case of the head-on collision of two BHs with 
unequal masses. We adopt therefore a coordinate system 
(r, 6, cp) adapted to the horizon T~L + so that r = const, char- 
acterizes sections St and we can write Sj = MDir, with 
A/" 2 = D i rD i r (i.e., M~ 2 = 7 rr ). Then, taking advan- 
tage of the axisymmetry, we adopt on St the preferred coor- 
dinated system (6, ip) discussed in Sec. IHI Al and consider the 
Cartesian-like coordinates constructed from (r, 6,<p) by stan- 
dard spherical coordinates relations: x = r sin6* coscp, y = 
r sm6 simp, z = r cosO. In these coordinates we have = 
A/(sin0 cos<^, smO simp, cos6*). Assuming the z-axis to be 
adapted to the axisymmetry, we choose £ l in Eq. d37l i as 
(£z) = M- 1 (d z ) i , so that (£, z ) l Si = cos9. Inserting ex- 
pression (|2"T1) in Eqs. d37b and d38l we obtain 



KT{t) = K°*U{t) 



(39) 



= "iF / ( cos ^) (f]^(t')Y efi (9)dA dA, 



with 



Af t (t)= / dt'h(t')+Ml 



(40) 



being the coefficients of a multipolar expansion of Eq. ( 1381 . 
Inserting the form dA = R^smOdOdip of the area element on 
St and performing the angular integration we finally find 



16tt 



(41) 
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with 



D 



(o) 

1,0 



V(2*+l)(2*-l) 



(42) 



As for the definition of A' h in the RT case, Eq. (HTt is 
quadratic in the (geometric) mass multipoles, i.e., the spher- 
ical harmonic components of the intrinsic curvature Ricci 
scalar curvature 2 R, although it involves a time integration, 
absent in (0. Also, it is an odd function under reflection with 
respect to z = const, planes and it involves only products 
of odd and even multipoles, precisely one of the criteria for 
the construction of K e g leading to the ansatz in Eq. (Q]0. In 
essence, expression (|4TT > for fulfills the two basic require- 
ments for the curvature parameter K c g with the added value 



that it is fully general and, in contrast with (O, no phenomeno- 
logical parameters need to be fitted. An additional and crucial 
feature is that terms involving I = 0, 1 are absent, due to the 
vanishing of I\ as discussed after Eq. (|2TT i. 

The quantity at the horizon H + is to be correlated with 
the component (dPf / dt) (t) of the flux of Bondi linear mo- 
mentum at J? + , which is useful to express in its multipolar 
expansion. First, we decompose $4 into its multipoles 



E 

£>2,m<\£\ 



(43) 



where ^■^T e ' m (6, <p) are the spin-weighted s = —2 spheri- 
cal harmonics. The explicit expression for the component of 
{dPf 1 / dt){t) along the z-axis (e.g. Ref. ll67tl ) then becomes 



dPi , , ,. r 2 
—r-(t) = lim — 



E ^ 

t>2,m<\l\ 



,(-2) rrl-1, 



(-2) Kfi+1, 



DT ' AT' 



(44) 



with 



dt'm 



being the corresponding multipolar components of the news 
functions introduced in ([34}, with the coefficients c\ ^ and 
D i7m given by 



C. 



(-2) 



2m 



1) 



D 



(-2) = 1 / (I - 2)(l + 2)(l-m)(l + m) 

£,m — ff\ 



(2£ - 1)(2£ + 1) 



(47) 



(45) 

The axisymmetric reduction of expression (144-b is obtained by 
setting m = in the expressions above. Note that ^4 is purely 
real in this cas^J The resulting coefficients are therefore 



(-2) 

to 



D 



(-2) 



{£-2){£ + 2) 



ro y {21- 1)(2£+ 1) 

(46) and we can write Eq. (1441 1 as 
I 



(48) 
(49) 



dPl 
dt 



(50) 



Expression ( l50b has obvious similarities with Eq. (HTt for 
K^ s (t). First, the (real) modes play a role analogous to 
those of the mass multipoles Ig. The common geometric na- 
ture of the underlying quantities ^4 and 2 R as curvatures, in 
particular, their dimensions as second derivatives of the met- 
ric, is indeed at the heart of the definition of the geometric 
multipoles I/s by Eqs. ( f20b and ( TZTT ) as the correct analogues 

H 

of Wf. Second, modes £ = 0, 1 are absent in both expres- 



7 Note that expression 44 1 1 cannot be factorized as a product of even and odd 
functions, as proposed in (TJ. 

8 The function JV rtc (Q) in Eq. 438t does not introduce £ = 1 modes either. 

9 For instance, the GW cross polarization h x vanishes: $4 = — H+ + ih x . 



I 

sions. This is nontrivial since the reasons underlying each 
case are different: the s = 2 spin weight of $4 in Eq. (l50l l 
and the vanishing of 1\ in d4TI ), respectively. This is a cru- 
cial feature for it directly impacts the determination of the 
mode dominating the dynamical behavior and, therefore, sin- 
gles out the Ricci scalar 2 R as a preferred quantity to be mon- 
itored instead of any other (spin-weighted s 7^ 2) function 
that could measure in some way the deformations of the hori- 
zon (for instance, the mean curvature). Besides the similari- 
ties, there are also differences between expressions d50b and 

(|4TT> . First, the coefficients Dfl in (|42| and (gSJl differ due 
to the different spin weight of 2 R and ^4. Therefore, the cor- 
relation between (dPf / dt)(t) and Kf s encodes information 
about the relative weight of the different couplings. Second, 
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the lower time-integration bound (t — s- — oo) is well-defined 
for (dPf /dt)(t), whereas (t) can be measured only after 
the formation of the common horizon. Finally, due to the ab- 
sence in the general case of a preferred coordinate system on 
St and their associated spherical harmonics, there is no nat- 
ural multipolar expression for K° s in the nonaxisymmetric 
case and one must resort to the full expression d3Tb . 

C. Correlation between the screens 

The effective-curvature vector introduced in the previ- 
ous section can now be used as a probe of the degree of cor- 
relation between the geometry at the horizon and the geom- 
etry far from the BH. More specifically, we aim at assessing 
the correlation between hi nn (t) — K^ s (t) at the horizon and 
hout{t) = {dPf / dt)(t) at large distances, considering these 
two quantities as time series. As discussed in Sec. [II] the use 
of a common time variable t for functions h nm and h out as- 
sumes a (gauge) mapping (cf. footnote 2 in Sec. [n} between 
the advanced v and retarded u times, parameterizing H + and 
J? + , respectively. The 3 + 1 slicing by hypersurfaces 
provides such a mapping, though an intrinsic time-stretching 
ambiguity between the signals at the two screens is present, 
due to the gauge nature of the slicing. This will be discussed 
in more detail later in this section. 

To quantify the similarities in the time series we employ 
the correlation function between time series h\(i) and h 2 (t), 
C(h\, h 2 ; t), defined as 

/oo 
hx(t + T)h 2 (t)dt. (51) 
-oo 

The structure of C(h\, h 2 ; t) encodes a quantitative compari- 
son between the two time series as a function of the time shift 
r (referred to as "lag") between them. This correlation func- 
tion encodes the frequency components held in common be- 
tween hi and h 2 and provides crucial information about their 
relative phases. Because the time series are intrinsically dif- 
ferent by a time lag, we measure the correlation between hi 
and h 2 as 

M(h u h 2 ) = max ( C Ahhhlll _) . (52) 

T \{C(hi,hi;0)C(h 2 ,h 2 ;0)}i J 

This number is confined between and 1 (where 1 indicates 
perfect correlation, and no correlation at all) and provides 
the maximum matching between the time series hi and h 2 
obtained by shifting one with respect to the other in time, and 
then normalized in frequency space. Besides providing a mea- 
sure of the correlation, expression Eq. d52l also gives a quanti- 
tative estimate of the coordinate time delay r max between the 
two signals. 

Note that one should not expect a perfect match between 
(dPf/dt){t) at and Kf (t) at H + , even if the latter re- 
sults to be a good estimator of the former. Indeed, (nonlinear) 
gravitational dynamics in the bulk spacetime affect and distort 
the possible relation between both quantities. Furthermore, 



given the related but different nature of 2 R and ^4 it is not 
obvious that a correlation should be found at all. 

In order to assess the validity of the approach, we con- 
struct Kl s {t) and (dPf / dt)(t) from the numerical simula- 
tions described in Sec. IIV Al Note that because Ii vanishes 
identically, the contributions N$Ni and AfiAf 2 are absent in 
the expression for K° s (t). Furthermore, since higher-order 
multipoles If become increasingly difficult to calculate, we 
truncate expression fiTt at I = 6; in our case, this has little 
influence on the overall results as we will show that the lowest 
(i.e., £ = 2,3) modes are by large the dominant ones. 

The values for (dPf /dt) and K° s as functions of the time 
t, and corresponding to the numerical simulations described in 
the previous Sec. IIV Al are presented in Fig. [12] The signals 
have been normalized with respect to their maximum value [a 
global rescaling does not affect the cross-correlations prop- 
erties of two functions hi(t) and h 2 (t), cf. Appendix [All . 
In Fig. [12] (a), the quantity K° is shown from the time 
t c pa 49.2M of first appearance of the common horizon (red 
dashed line). After time £ max w 120M the error in the calcu- 
lation of the I( multipoles becomes comparable with the value 
of the multipoles, spoiling the evaluation of the integrals in 
d40l ). Hence, we set to zero for t > £ max . Similarly, the 
flux of Bondi linear momentum {dPf / dt) as computed by an 
observer at 100 M from the origin, is split in a part before the 
appearance of the common AH (blue dotted line) and in one 
which is to be compared with (blue dashed line). In pan- 
els (b) and (c) of Fig. [12] we show instead K ° s and (dP® /dt) 
separately and in different time intervals for a better emphasis 
of the similarities. 

Some interesting remarks on Fig. [12] can be made already 
at a qualitative level. In particular, it is clear that suc- 
ceeds in tracking key features of (dP® /dt). This is apparent 
in the relative magnitude of the three first positive peaks in 
the two signals and the qualitative agreement is maintained in 
time. As expected, some specific features of (dPf /dt) are not 
faithfully captured in , such as the magnitude of the nega- 
tive peak around t ps 148M relative to the neighboring peaks. 
However given the heuristic character of K ° ff and the fact that 
its geometric definition does not leave room for any tuning, 
the overall qualitative agreement with {dP® /dt) at al- 
ready represents a remarkable result, shedding light on the 
near-horizon dynamics. This agreement between (dPf/di) 
and A'° ff is indeed the main result of this section and the ulti- 
mate justification for the introduction of K° s . It is also worth 
stressing that attempts employing other quantities (e.g. a blind 
application of the methods used for RT spacetimes) would not 
lead to such a clear matching. 

From a quantitative point of view, the correlation analysis 
for the time intervals shown in Fig. Q~2] (b) and (c) indicates 
that the two signals yield a typical correlation Ai « 0.93 and 
a time lag r = 97 M (we recall that the observer is at 100 M 
and that the common AH has the size of a couple of M). How- 
ever, as one tries to extend the analysis to the very first time 
of the formation of the apparent horizon, the correlation drops 
significantly. The reason for this drop is related to the stretch- 
ing of the time coordinate between the two screens. In addi- 
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FIG. 12. Effective curvature Kf 1 defined at the horizon via Eq. J37l > (red dashed curve) and flux of momentum (dPf /dt) evaluated at an 
approximation of with Eq. d33 b (blue dotted and solid curves). These quantities encode, respectively, the information of the common 
horizon deformation and the flux of momentum carried away by GWs in the head-on collision of BHs with mass ratio 1/2. Note, at panel (a), 
the qualitative agreement between those curves, which allows us to distinguish the momentum radiated before (blue dotted curve) and after 
(blue solid curve) the merge. Panel (b) and (c) compares the same quantities for latter times, where one can still see the good agreement. 



tion to the obvious time delay between the (dPf / dt){t) and 
K° s (t) due to the finite (coordinate) speed of light, in fact, 
the dependence of the two signals in coordinate time t is not 
the same and is stretched between the two screens. This effect 
is the result of the in-built gauge mapping between sections of 
J? + and the horizon H + defined by the spacetime slicing, but 
also of the physical blueshift (redshift) of signals at the inner 
(outer) boundaries in the BH spacetime. 

Although approaches to disentangle the physical and gauge 
contributions can be derived, for instance by introducing 
proper times of suitably defined observers, this goes beyond 
the scope of this paper. Rather, we opt here for a more 
straightforward approach in which only comparisons based 
on sequences of (absolute values and signs of the) maxima 
and minima in the signals h- nm (t) and h out (t) are considered 
significant, since the relative shape of h- mn (t) and h out (t) can 
be subject to a time reparametrization. This association is pos- 
sible when the quantities which are compared are scalars, so 
that the values of maxima and minima are well-defined and 
independent of coordinates. This is possible in the case of 
axisymmetry as it gives a privileged direction along which 
to contract the effective-curvature vector Kf . In a more 
generic configuration one would need to build an appropriate 
frame to produce scalars by contraction with tensorial quanti- 
ties. Once the correspondence between maxima and minima 
in the two signals hi nn (t) and h out (t) is established, a map- 
ping t out = ^out(imn) can be easily constructed. With this 
matching, the calculation of the correlation parameter gives 
typically values A4 ^ 0.9 for any chosen time interval 0- 
More information about the mapping between t out and t- ln is 



found in AppendixlBl 

1. A critical assessment of the correlation 

Of course, it is reasonable to question if finding such a high 
correlation is just a bias in our methodology. Certainly, our 
strategy of identification of maxima and minima in the sig- 
nals enhances the correlations when time-stretching issues are 
involved. However, it does not guarantee by itself the high 
(positive) values found for Ai. More specifically, once a 
first couple of maxima (or minima) are identified in the two 
signals, the subsequent couples of maxima and minima con- 
structed from the data in each signal are automatically fixed. 
As a consequence, high correlations are possible only if the 
sequence of signs in the extrema of the two signals is exactly 
the same. In addition to consideration above, one may also 
argue that the high correlation found is just the result of the 
very rapid decay of the signals, which makes the first couple 
of maxima and minima play a dominant role in the estimate, 
possibly shadowing the role of the smaller peaks appearing 
at later times. To address this point and weight equally all 
parts of the signals, we model them as exponentially decay- 
ing oscillating functions, i.e., hf nn (t) = e Kinnt /i; nn (i) and 
^out(^) = eKout *^out(i)- This is applied to the signals without 
the time correction provided by the map t out = t ou t (ti nn ) , 
finding 

Mn inn = 0.179 ± 0.005 , Mn out = 0.181 ± 0.006 , 

(53) 



Interesting information can also be gained by studying in more detail the 
properties of the mapping t ou t = tout(iinn)- More specifically, we have 
found that the derivative dt ou t /dti nn is not constant and starts as being 



larger than unity (indicating that initially the coordinate time at runs 
faster than the time at W"""), but then oscillates around unity at late times. 
This is consistent with the fact that as stationarity is approached, the evo- 
lution vector t a adapts to the timelike Killing vector. 
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FIG. 13. Quantities on H + (red dashed line) and ^ + (blue solid 
line) as shown in Fig.[l2]but without an overall exponential decay in 
time. The close similarity in the signals is confirmed by the very large 
correlation which is A4 — 0.87, obtained without a time mapping. 



through a least-square fitting. The resulting functions hf nn (t) 
and h% ut (t) once the exponential decay has been subtracted 
are shown in Fig. [13] Once again it is apparent that the 
two time series are very similar and indeed the matching 
computed even without introducing any time mapping is 
M(h* nn , h* ut ) = 0.87 and thus remarkably high. 

The main reason behind the good correlation found also for 
the undamped signals is that the next-to-leading-order term, 
i.e., term A/3A/4 and the corresponding i = 3 and I = 4 cou- 
pling in Eq. ( T50t . are much smaller than the leading-order term 
A2A/3. Indeed, we have found that it is possible to express to 
a very good approximation K° s ~ M2M3 and (dPf /dt) ~ 
A2A/3. This is confirmed by the corresponding power spectra, 
which are shown in Fig. [14] and are dominated in both cases 
by two frequencies: Clf n = 0.22 ± 0.04, 2 nn = 0.98 ± 0.05 
for the signal h inn (t), and ft? ut = 0.22 ± 0.04, Qf lt = 
0.97 ± 0.04 for the signal h ont (t). 

These frequencies are closely related to the quasinormal 
modes of the merged BH, interfering to lead to a beating sig- 
nal. To see this, we model each function fifi and Me as an 



FIG. 14. Normalized power spectrum for (fi) (red dashed line) 
and {dP z / dt)(Q) (blue solid line) as measured at H + and 
respectively (cf. Figure 1 1 3b - Both spectra are dominated by two 
frequencies: f)f n = 0.22 ± 0.04 and nif n = 0.98 ± 0.05 (red 
dashed line) and f)? ut = 0.22 ± 0.04 and * = 0.97 ± 0.04 (blue 
solid line) which are linear combinations of the quasinormal ringing 
modes of the merged BH (cf. Table[l](. 



at JP + , consistent with the "beating" behavior shown by the 
power spectra in Fig. [14] Similarly, the decay time scales are 
then given by 



= 3 ! 



K e=3- 



(56) 



These frequencies and time scales match very well the 
real (uif') and imaginary (uj l e ) parts of the fundamental 
(n = 0) quasi-normal-modes (QNM) eigenfrequencies of a 
Schwarzschild BH l68ll . A detailed comparison is presented 
in Table [IJ whose first six columns report the properties of 
the signals hi nn (t) and h out (t) in their constituent frequen- 
cies 3 and 0^1.3 3 defined in Eqs. (f54b— (T55l). and com- 
pare them with the corresponding real parts of the eigenfre- 
quencies of a Schwarzschild BH, w|L 2 3 - The close match in 
the oscillatory part is accompanied also by a very good cor- 
respondence in the decaying part of the signal. Defining, in 
fact, the overall decay time in terms of the imaginary parts of 



exponentially damped sinusoid, i.e., Me ~ e _K « t sm(flf t + the QNM eigenfrequencies, i.e., as Kdocay = 



,e=3 



■r 



and Me — e _K « *sin(^ t + yA). Then, under the ap- 
proximation K° s — M2M3 and {dPf /dt) — M2M3 it follows 



1 If — ' 



n 2 nn - n 1 / 1 



^fc3 



at T~L + , whereas 



QOUt (~)OUt 

a M _ "2 "1 W 

lL e=2 — o ' "£=3 



fl 2 nn + ^i n 



Qout + fi out 



(54) 



, (55) 



is easy to realize from the last three columns in Table [I] that 
this decay time is indeed very close to the one associated to 
the signal at the two screens [cf. Equations $55[ and (l56l)l. 

This role of QNMs is not entirely surprising for a measure 
at J? + , but it is far less obvious to see it imprinted also for a 
quantity measured at H + . This indicates that the bulk space- 
time dynamics responsible for the recoil physics is a relatively 
mild one, so that a QNM ringdown behavior dominates the dy- 
namics of the deformed single AH and imprints the properties 
of the radiated linear momentum. It is interesting that a purely 
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(quasi-)local study of the AH geometric properties permits us 
to read the behavior of quantities which are intrinsically de- 
fined at infinity, thus confirming the main thesis in Ref. QJJ] . 



2. Antikicks and the Slowness Parameter 

As a concluding remark for this section, we make use of our 
results, and, in particular, on the spectral and decaying prop- 
erties of our measures on the screens, to make contact with 
the analysis carried out in l32tl . More specifically, we can 
define a characteristic decay time r = (27r)/K mn / out and an 



inn/out 



from which 



oscillation-characteristic time T = 2n/n. 
to build our equivalent of the "slowness parameter" P = TJt 
introduced in l32tl . The specific case discussed above then 
yields r ~ 34.9M, T ~ 6AM and thus P ~ 0.18. As de- 
tailed in Ref. l32ll . small antikicks should happen when the 
two timescales are comparable, thus corresponding to an os- 
cillation which is over-damped. This expectation is indeed 
confirmed by the recoil velocity shown in Fig. [TO] where the 
relative antikick is about ~ 30% and thus compatible with 
the slowness parameter that we have associated to our process 
[see also the discussion below on the application of Eqs. (T57T > 
and ( f58l > in Fig. [15). This qualitative agreement with the phe- 
nomenological approach discussed in Ref. l32ll is very natu- 
ral. While we here concentrate on modeling the local curva- 
ture properties at the horizon, Ref. l32tl concentrates on the 
spectral features of the signal at large distances. Since we 
have demonstrated that the two are highly correlated, it does 
not come as a surprise that the two approaches are compati- 
ble. Looking at the local horizon's properties has however the 
added value that it provides a precise framework in which to 
predict not only the strength of the antikick, but also its di- 
rectionality. Furthermore, such an approach permits an inter- 
pretation of BH dynamics in terms of viscous hydrodynamics, 
as we will discuss in detail in paper II. In particular, we shall 
show there that the horizon-viscous analogy naturally leads to 
a geometric prescription for an (instantaneous) slowness pa- 
rameter P, in terms of timescales r and T respectively related 
to bulk and shear viscosities. 

The logic developed above for the calculation of the slow- 
ness parameter can be brought a step further by assuming that 
the final BH produced by the merger of a binary system in 
quasicircular orbit can be described at the lowest order by an 
oscillation and decay times 



2n 



T 



2 71 



UJ 



J e=3 



U f=2 + U f= 



to which corresponds a slowness parameter defined as 

p _ X. = ^4=2 + ^1=3 



(57) 



(58) 



1=2 



Using the semianalytic expressions derived for estimating the 
spin of the final BH, e.g. I69t - l7lll . it is possible to predict 
the values of r and T for any binary whose initial spins and 
masses are known, and thus predict qualitatively through P 
the strength of the antikick which will be produced in any of 



these configurations. We have tested this idea by consider- 
ing the data presented in Ref. lf20ll both for the kick/antikick 
velocities and for the final spin of the merged BH. This con- 
jecture about the predictability of the antikick in terms of 
the slowness parameter is indeed supported by the example 
data collected in Fig. Q3] More specifically, the left panel in 
Fig. [TBI shows the correlation between the slowness parame- 
ter P = T/t as computed in Eq. ( l58b and the dimensionless 
spin of a BH afi n = J^a/M^ produced, for example, in the 
merger of a binary system (expressions to estimate the QNM 
eigenfrequencies for rotating BHs can be found in a number 
of works which are collected in the review (6§H ). The middle 
panel shows instead the good correlation between the relative 
antikick velocity Av/v &11 = (w max - v fin )Mn = Vk/vn„ and 
the dimensionless final spin as computed from the data taken 
from Ref. ll2bTl (indicated with error bars are the estimated nu- 
merical errors). Finally, the right panel combines the first two 
and shows the searched correlation between the antikick ve- 
locity and the slowness parameter. It also shows with a solid 
line an exponential fit, which suggests a vanishing antikick for 
a slowness parameter P ~ 1. All in all, this figure confirms 
also for the case of binaries in quasicircular orbits the sug- 
gestion 13211 that the smaller the slowness parameter P gets, 
the larger is the expected value of the antikick. Large anti- 
kicks should then be expected for F< 1 ll32Tl . Furthermore, 
it highlights that it is indeed possible to predict qualitatively 
the antikick merely on the basis of the initial properties of the 
BHs when the binary is still widely separated. 



V. CONCLUSIONS 

We have demonstrated that qualitative aspects of the post- 
merger recoil dynamics at infinity can be understood in terms 
of the evolution of the geometry of the common horizon of 
the resulting black hole. This extends to binary black-hole 
spacetimes the conclusions presented in Ref. [Q]] based on 
Robinson-Trautman spacetimes. More importantly, we have 
shown that suitably built quantities defined on inner and outer 
world tubes (represented either by dynamical horizons or by 
timelike boundaries) can act as test screens responding to the 
spacetime geometry in the bulk, thus opening the way to a 
cross-correlation approach to probe the dynamics of space- 
time. 

The extension presented here is nontrivial and it involves 
the construction of a phenomenological vector Kf s (t) from 
the Ricci curvature scalar 2 R on the dynamical-horizon sec- 
tions, which then captures the global properties of the flux 
of Bondi linear momentum (dP® / dt) (t) at infinity, namely, 
(proportional to) the acceleration of the BH. At the same time, 
the proposed approach involves the development of a cross- 
correlation methodology which is able to compensate for the 
in-built gauge character of the time evolution on the two sur- 
faces. A proper mapping between the times on the two sur- 
faces is needed and its gauge nature highlights that the phys- 
ical information encoded in the surface quantities is not in 
its local (arbitrary) time dependence, but rather in the global 
structure of successive maxima and minima. 
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TABLE I. The first six columns offer a comparison between the properties of the signals hi nn (t) and h ou t(t) in their constituent frequencies 
f2^C 2 ,3 defined in Eqs. J541> — J55b. with the corresponding real parts of the eigenfrequencies of a Schwarzschild BH, w|L 2 3 . The last three 
columns show instead a comparison between the damping times /thin.out defined in Eq. d56t . with the corresponding decay time ftdecay 
computed from the imaginary parts of the eigenfrequencies. In all cases, the close match is remarkable and not at all obvious for quantities 
measured at H + . 
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-^^decay 


0.38 ±0.04 


0.37 ±0.04 
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0.60 ± 0.04 
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0.181 ±0.006 
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FIG. 15. Predictability of the antikick in terms of the slowness parameter and the antikick. The left panel shows the correlation between 
the slowness parameter P = T/t as computed in Eq. {58} and the dimensionless spin of a BH an n . The mid panel shows instead the good 
correlation between the relative antikick velocity Av/va n — Uk/ufln and the dimensionless final spin, using the data taken from Ref. (20ll 
(indicated with error bars are the estimated numerical errors). Finally, the right panel combines the first two and shows the correlation between 
the antikick velocity and the slowness parameter. In particular, the good exponential fitting shows the consistency with a vanishing antikick 
for a slowness parameter P ~ 1, as discussed in (32ll . 



By analyzing Robinson-Trautman spacetimes, Ref. |l|] pro- 
posed that when a single horizon is formed during the merger 
of two BHs, the observed decelerations/accelerations of the 
newly formed BH can be understood in terms of the dissipa- 
tion of an anisotropic distribution of the Ricci scalar curva- 
ture on the horizon. The results presented here confirm this 
picture, although through quantities which are suited to BH 
spacetimes. Being computed on the horizon, these quantities 
reflect the properties of the BH and, in particular, its exponen- 
tially damped ringing. The interplay between oscillation and 
decay timescales associated with this process, which are in- 
evitably imprinted in our geometric variables, explain the late 
qualitative features of the recoil dynamics, in particular, the 
antikick, in natural connection with the approach discussed in 
Ref. 13211 . where the antikick is explained in terms of the spec- 
tral features of the signal at large distances. Because we have 
shown that the latter is closely correlated with the signal at 
the horizon, we can adopt the same slowness parameter intro- 
duced in Ref. l32ll to predict qualitatively the magnitude of the 
antikick from the merger of BH binaries with spin aligned to 



the orbital angular momentum, finding a very good agreement 
with the numerical data. 

As a final remark we note that looking at the horizon's prop- 
erties has the added value that it provides a precise frame- 
work in which to predict not only the strength of the antikick, 
but also its directionality. Furthermore, as we discuss in de- 
tail in paper II, our geometric (cross-correlation) framework 
presents a number of close connections with (and potential 
implications on) the literature developing around the use of 
horizons to study the dynamics of BHs, as well as with the 
interpretations of such dynamics in terms of a viscous hydro- 
dynamics analogy. Much of the machinery developed using 
dynamical trapping horizons as inner screens can be extended 
also when a common horizon is not formed (as in the calcula- 
tions reported in Ref. |31j). While in such cases the identifica- 
tion of an appropriate hypersurface for the inner screen can be 
considerably more difficult, once this is found its geometrical 
properties can be used along the lines of the cross-correlation 
approach discussed here for dynamical horizons. 
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We note that, according to expression ( lAU . the numerator in 
the second line of Eq. ( IA5I ) is the inverse Fourier transform 
of the cross-spectrum function C{h\, h 2 )(f). Therefore the 
numerator in the expression for JA(hi, h 2 ) is just the corre- 
lation function C(hi, h 2 )(r). Regarding the denominator, we 
use Parseval's identity 



/ Mf)\ 2 df= / \hi(t)\ 2 dt, (A6) 

J —oo J —oo 



Appendix A: Correlation and matching of time series 



and the expression for the autocorrelation of functions hi(t) 



The correlation function C(hi,h 2 ;r) introduced in 
Sec. HVCI provides information in the time domain about the 
comparison of temporal series hi(t) and h 2 (t). Its Fourier 
transform defines the cross-spectrumC(hi,h 2 ; f) of hi and 
h 2 , providing the corresponding analysis in the frequency do- 
main. It has the form 



C(h 1 ,h 2 ;f) = h 1 (f)h* 2 (f), 
where Fourier transform conventions are 



(Al) 



hi(f) 



hi{t)e i2nft dt , hi(t) 



hi(f)e- a " ft df . 

(A2) 

Choosing a measure (S n (\f\)) 1 df, a natural scalar product 
between functions h\(t) and h 2 (t) (or h\(f) and h 2 (f)) is 
introduced as 



(hi, hi 



hi(mu) 

Sn{\f\) 



df. 



(A3) 



In GW data analysis S n (\f\), the noise power-spectral den- 
sity, is associated with the spectral sensitivity of the instru- 
ment. In our case we have no a priori knowledge about 
S n (\f\), and we choose S n (f) = 1. The scalar product (•, •) 
introduces the natural projection between hi(f) and h 2 (f). 
Their normalized scalar product defines the overlap 



0[hi,h 2 



(hiM) 



(hi, hi)(h 2 ,h 2 ) 



(A4) 



Fixing one of the functions, say hi(t), we can consider its 
overlap with the function resulting by shifting h 2 (t) in time 
by a time lag r, i.e., h 2 (t + r). In the frequency do- 
main this amounts to calculate the overlap between hi(t) and 
h 2 (f,r) = h 2 (f)e~ t2 ' K * T . Maximizing over r provides the 
best match estimator 



A4(hi, h 2 ) = max{C[/ii, h 2 (r)]} = max 



(hi,h 2 (r)) 



max ■ 



\J(hi,hi)(h 2 ,h 2 ) 



(j-oeMfiMJZMfWdf 



(A5) 



\hi(t)\ 2 dt = C(hi,hi;T = Q) 



(A7) 



We then recover expression (l52t for M(hi, h 2 ) in terms of 
the correlation function C( h i, h 2 ;r). 



Appendix B: Mapping time series on the screens 



As discussed in Sec. IIV CI a built-in gauge mapping be- 
tween sections of J? + and the horizon H + defined by the 
spacetime slicing leads to a stretching of the time coordinate 
between the two screens. 

A comparison based on sequences of maxima and minima 
in the signals h- mn (t) and h out (t) allows us to construct the 
mapping t ou t(t- mn ), which is here depicted in FigfTSh. Also 
shown in FigfTSb is the derivative of this interpolation, crucial 
to assess the relative rate of the considered coordinate times. 
In particular, it addresses the behavior discussed in footnote 

Coli. 

Finally, FigJT6b presents the correlation number M. be- 
tween the two signals as a function of time intervals A. The 
construction of the intervals A is based on the sequences of 
maxima and minima identified in the two signals. In this 
way, we fix the final time in both series as tf n "j al = 96. 8M 



final 
out 

t final _ t initial ) starting f rom ^initial = 49.2M and t™^ 1 



and = 194.471/ and then we establish windows A 

^initial 
inn 



140.4M. For the red curve, the correlation is evaluated with- 
out using the mapping t out = t ut(£mn) t° correct the stretch- 
ing, while the black curve takes the effect into account. This 
latter figure shows that the correction through the mapping 
tout = tout(tinn) is crucial to disentangle coordinate from 
real effects at early times. 
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FIG. 16. Mapping between the coordinate time t[ n and t out measured at the inner/outer screens. The first panel shows the interpolation of the 
function i out = iout(£in) constructed from the comparison between the sequence of maxima and minima in the signals hi nn (t) and h ou t(t). 
The middle panel depicts the interpolation's derivative dt a ut/dt\ n . In particular, it shows that initially the coordinate time at runs faster 
than the time at H + and than oscillates around unity at late times. This behavior is consistent with the approach to stationarity. Finally, the 
right panel presents the correlation number Ai as a function of time intervals A. The red curve shows the correlation without the using the 
mapping i ou t = tout (tin) to correct the time stretching between the signals on the two screens, whereas the black curve takes the effect into 
account and gives M ^ 0.9. This correction is crucial to disentangle coordinate from real effects at early times. 
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